Question: Does a continuous wetland to upland gradient exist? Or, can a catena-type model be derived using wetlands and uplands in a data-driven approach?

If the ends of the wetland/upland classification are accurate, does the intermediate range correspond with hypothesized soil formation characteristics according to a hydropedological framework?

Hydropedological units in the Hubbard Brook Experimental Forest include:

hydropedology hillslope model
hydropedology hillslope model

Using Hubbard Brook as an example



hbdem <- rast("UplandWetlandGradient/data/dem1m.tif")
plot(hbdem)

hb <- vect("UplandWetlandGradient/data/hbef_wsheds/hbef_wsheds.shp")

lines(hb)

nwi <- vect("UplandWetlandGradient/data/HU8_01070001_Watershed/HU8_01070001_Wetlands.shp") |>
    project("EPSG:26919")
hbnwi <- nwi |>
    crop(ext(hbdem))  #|> dplyr::filter(WETLAND_TY != 'Riverine')
names(hbnwi)
## [1] "ATTRIBUTE"  "WETLAND_TY" "ACRES"      "Shape_Leng" "Shape_Area"
plot(hbdem)
plot(hbnwi, "WETLAND_TY", type = "classes", add = T)

hbgeo <- vect("UplandWetlandGradient/data/hbef_bedrock/hbef_bedrock.shp")


# 1 meter LiDAR-derived, Filtered Hydro-enforced Digital
# Elevation Model (DEM)
hb1mlpns <- rast("UplandWetlandGradient/data/hydem1mlpns.tif")


plot(hb1mlpns)
plot(hbgeo, type = "classes", add = T)

Merge hydrography and NWI

hbhydrograph <- vect("UplandWetlandGradient/data/hbef_hydro/hbef_hydro.shp") |>
    terra::buffer(0.5)
plot(hbhydrograph, type = "classes")

hbnwi_hydrog <- hbnwi |>
    terra::aggregate("WETLAND_TY") |>
    terra::union(y = hbhydrograph)
plot(hbnwi_hydrog, type = "classes")

Get training data - wetland and upland points

set.seed(11)
randupl <- spatSample(hbdem, 2000, as.points = TRUE, xy = TRUE)

hbnwi_buff <- buffer(hbnwi_hydrog, 10)

hbupl_pts <- terra::mask(randupl, hbnwi_buff, inverse = T) |>
    mutate(class = "UPL") |>
    select(-dem1m)


plot(hbupl_pts, "class", type = "classes")
lines(hbnwi_buff)

hbwet_pts <- spatSample(hbnwi_hydrog, size = 500) |>
    mutate(class = "WET") |>
    select(class)

plot(hbwet_pts, "class", type = "classes", col = "blue", cex = 0.2)
plot(hbupl_pts, "class", type = "classes", col = "red", cex = 0.2,
    add = T)

hbpts_all <- rbind(hbupl_pts, hbwet_pts)
hbpts_ext <- terra::extract(hbdem, hbpts_all, bind = T)
hbpts <- hbpts_ext |>
    dplyr::filter(!is.na(dem1m))
writeVector(hbpts, "UplandWetlandGradient/data/derived_data/hbpts_hydrog_raw.gpkg",
    overwrite = TRUE)

plot(hbdem)
plot(hbhydrograph, type = "classes", add = T)
plot(hbpts, "class", type = "classes", cex = 0.3, add = T)

Make terrain metrics

hbslp_3 <- SlpAsp(hbdem, w = c(3, 3), metrics = "slope", filename = "UplandWetlandGradient/data/derived_data/hbslp_3.tif",
    overwrite = T)

hbslp_27 <- SlpAsp(hbdem, w = c(27, 27), metrics = "slope", filename = "UplandWetlandGradient/data/derived_data/hbslp_27.tif",
    overwrite = T)

hbslp_81 <- SlpAsp(hbdem, w = c(81, 81), metrics = "slope", filename = "UplandWetlandGradient/data/derived_data/hbslp_81.tif",
    overwrite = T)

hbtpi_3 <- TPI(hbdem, w = c(3, 3), shape = "rectangle", stand = "none",
    na.rm = TRUE, filename = "UplandWetlandGradient/data/derived_data/hbtpi_3.tif",
    overwrite = T)

hbtpi_27 <- TPI(hbdem, w = c(27, 27), shape = "rectangle", stand = "none",
    na.rm = TRUE, filename = "UplandWetlandGradient/data/derived_data/hbtpi_27.tif",
    overwrite = T)

hbtpi_81 <- TPI(hbdem, w = c(81, 81), shape = "rectangle", stand = "none",
    na.rm = TRUE, filename = "UplandWetlandGradient/data/derived_data/hbtpi_81.tif",
    overwrite = T)

hbcurv_3 <- Qfit(hbdem, w = c(3, 3), metrics = c("meanc", "profc",
    "planc"), unit = "degrees", na.rm = TRUE, filename = "UplandWetlandGradient/data/derived_data/hbcurv_3.tif",
    overwrite = T)

hbcurv_27 <- Qfit(hbdem, w = c(27, 27), metrics = c("meanc",
    "profc", "planc"), unit = "degrees", na.rm = TRUE, filename = "UplandWetlandGradient/data/derived_data/hbcurv_27.tif",
    overwrite = T)

hbcurv_81 <- Qfit(hbdem, w = c(81, 81), metrics = c("meanc",
    "profc", "planc"), unit = "degrees", na.rm = TRUE, filename = "UplandWetlandGradient/data/derived_data/hbcurv_81.tif",
    overwrite = T)

Slope at different scales

plot(hbslp_3)
plot(hbslp_27)
plot(hbslp_81)

TPI at different scales

plot(hbtpi_3)
plot(hbtpi_27)
plot(hbtpi_81)

Curvature at different scales

plot(hbcurv_3)
plot(hbcurv_27)
plot(hbcurv_81)

super_stack <- c(hbslp_3, hbslp_27, hbslp_81, hbcurv_3, hbcurv_27,
    hbcurv_81, hbtpi_3, hbtpi_27, hbtpi_81)

hbpts_extract <- hbpts |>
    terra::extract(x = super_stack, bind = TRUE, xy = TRUE) |>
    drop_na()
names(hbpts_extract) <- (c("class", "dem", "slp_3", "slp_27",
    "slp_81", "meancurv_3", "prof_curv_3", "plan_curv_3", "meancurv_27",
    "prof_curv_27", "plan_curv_27", "meancurv_81", "prof_curv_81",
    "plan_curv_81", "tpi_3", "tpi_27", "tpi_81", "x", "y"))
hbpts_extract$class <- as.factor(hbpts_extract$class)

writeVector(hbpts_extract, "UplandWetlandGradient/data/derived_data/hbpts_hydrog_extract.gpkg",
    overwrite = TRUE)
hist(hbpts_extract$slp_3)
hist(hbpts_extract$slp_27)
hist(hbpts_extract$slp_81)

Slope at 3m vs. 81m

plot(hbpts_extract$slp_3, hbpts_extract$slp_81)
abline(0, 1, col = "red")

TPI at 3m vs. 81m

plot(hbpts_extract$tpi_3, hbpts_extract$tpi_81)
abline(0, 1, col = "red")

Mean Curvature at 3m vs. 81m

plot(hbpts_extract$meancurv_3, hbpts_extract$meancurv_81)
abline(0, 1, col = "red")

Is resampling the same as setting the window size larger?

hbcurve_agg <- terra::aggregate(hbcurv_3, fact = 81, fun = "mean")

hbpts_extract_agg <- hbpts_extract |>
    terra::extract(x = hbcurve_agg, bind = T) |>
    drop_na()
names(hbpts_extract_agg)
##  [1] "class"        "dem"          "slp_3"        "slp_27"       "slp_81"      
##  [6] "meancurv_3"   "prof_curv_3"  "plan_curv_3"  "meancurv_27"  "prof_curv_27"
## [11] "plan_curv_27" "meancurv_81"  "prof_curv_81" "plan_curv_81" "tpi_3"       
## [16] "tpi_27"       "tpi_81"       "x"            "y"            "meanc"       
## [21] "profc"        "planc"
hbpts_extract_agg
##  class       : SpatVector 
##  geometry    : points 
##  dimensions  : 2246, 22  (geometries, attributes)
##  extent      : 274043.4, 283958.1, 4866042, 4871958  (xmin, xmax, ymin, ymax)
##  coord. ref. : NAD83 / UTM zone 19N (EPSG:26919) 
##  names       :  class   dem slp_3 slp_27 slp_81 meancurv_3 prof_curv_3
##  type        : <fact> <num> <num>  <num>  <num>      <num>       <num>
##  values      :    UPL 333.2 14.31  12.14  9.286   -0.02777    -0.04499
##                   UPL 745.8 15.09  8.697  12.81   -0.03448     0.00313
##                   UPL   584 14.21  15.18  15.11   0.001648    0.007331
##  plan_curv_3 meancurv_27 prof_curv_27 (and 12 more)
##        <num>       <num>        <num>              
##     -0.01056   0.0003643   -0.0001417              
##     -0.07208  -0.0001351     0.008627              
##    -0.004036    0.002487     0.003541
plot(hbcurve_agg)

Set up random forest model

library(randomForest)
## randomForest 4.7-1.1
## Type rfNews() to see new features/changes/bug fixes.
## 
## Attaching package: 'randomForest'
## The following object is masked from 'package:dplyr':
## 
##     combine
## The following object is masked from 'package:spatialEco':
## 
##     combine
set.seed(11)

hbpts_exdf <- as.data.frame(hbpts_extract)

# Validation Set
train.index <- as.vector(sample(c(1:nrow(hbpts_exdf)), 0.7 *
    nrow(hbpts_exdf), replace = F))

train <- hbpts_exdf[train.index, c("class", "dem", "slp_3", "slp_27",
    "slp_81", "meancurv_3", "prof_curv_3", "plan_curv_3", "meancurv_27",
    "prof_curv_27", "plan_curv_27", "meancurv_81", "prof_curv_81",
    "plan_curv_81", "tpi_3", "tpi_27", "tpi_81")]

test <- hbpts_exdf[-train.index, c("class", "dem", "slp_3", "slp_27",
    "slp_81", "meancurv_3", "prof_curv_3", "plan_curv_3", "meancurv_27",
    "prof_curv_27", "plan_curv_27", "meancurv_81", "prof_curv_81",
    "plan_curv_81", "tpi_3", "tpi_27", "tpi_81")]
library(caret)
## Loading required package: ggplot2
## 
## Attaching package: 'ggplot2'
## The following object is masked from 'package:randomForest':
## 
##     margin
## Loading required package: lattice
rf_model <- randomForest((class) ~ ., ntree = 500, na.action = na.omit,
    importance = TRUE, data = train)
plot(rf_model)

rf.test <- predict(rf_model, newdata = test, type = "response")
rf.train <- predict(rf_model, newdata = train)
caret::confusionMatrix(rf.train, train$class)
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction  UPL  WET
##        UPL 1264    0
##        WET    0  309
##                                      
##                Accuracy : 1          
##                  95% CI : (0.9977, 1)
##     No Information Rate : 0.8036     
##     P-Value [Acc > NIR] : < 2.2e-16  
##                                      
##                   Kappa : 1          
##                                      
##  Mcnemar's Test P-Value : NA         
##                                      
##             Sensitivity : 1.0000     
##             Specificity : 1.0000     
##          Pos Pred Value : 1.0000     
##          Neg Pred Value : 1.0000     
##              Prevalence : 0.8036     
##          Detection Rate : 0.8036     
##    Detection Prevalence : 0.8036     
##       Balanced Accuracy : 1.0000     
##                                      
##        'Positive' Class : UPL        
## 
caret::confusionMatrix(rf.test, test$class)
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction UPL WET
##        UPL 528  25
##        WET  14 108
##                                           
##                Accuracy : 0.9422          
##                  95% CI : (0.9219, 0.9586)
##     No Information Rate : 0.803           
##     P-Value [Acc > NIR] : <2e-16          
##                                           
##                   Kappa : 0.8115          
##                                           
##  Mcnemar's Test P-Value : 0.1093          
##                                           
##             Sensitivity : 0.9742          
##             Specificity : 0.8120          
##          Pos Pred Value : 0.9548          
##          Neg Pred Value : 0.8852          
##              Prevalence : 0.8030          
##          Detection Rate : 0.7822          
##    Detection Prevalence : 0.8193          
##       Balanced Accuracy : 0.8931          
##                                           
##        'Positive' Class : UPL             
## 
vip::vip(rf_model)

AUC/ROC

library(performance)
library(pROC)
## Type 'citation("pROC")' for a citation.
## 
## Attaching package: 'pROC'
## The following objects are masked from 'package:stats':
## 
##     cov, smooth, var
library(ROCR)
## 
## Attaching package: 'ROCR'
## The following object is masked from 'package:performance':
## 
##     performance

rf_prediction_test <- ROCR::prediction(as.numeric(rf.test), as.numeric(test$class))
roc_test <- ROCR::performance(rf_prediction_test, measure = "tpr",
    x.measure = "fpr")
auc_test <- ROCR::performance(rf_prediction_test, measure = "auc")

plot(roc_test)
abline(0, 1, col = "red", lty = 5)
print(auc_test@y.name)
## [1] "Area under the ROC curve"
print(auc_test@y.values)
## [[1]]
## [1] 0.8930999

Predict maps

pred_stack <- c(hbdem, super_stack)


names(pred_stack) <- (c("dem", "slp_3", "slp_27", "slp_81", "meancurv_3",
    "prof_curv_3", "plan_curv_3", "meancurv_27", "prof_curv_27",
    "plan_curv_27", "meancurv_81", "prof_curv_81", "plan_curv_81",
    "tpi_3", "tpi_27", "tpi_81"))
hbwip <- predict(pred_stack, rf_model, type = "prob", filename = "UplandWetlandGradient/data/derived_data/hbwip_hydrog.tif",
    overwrite = TRUE)

plot(hbwip$WET)
lines(hb)

plot(hbwip_noriv$WET)
lines(hb)

PDF/kernel density plot

hbwip_trainpred <- predict(rf_model, newdata = train, type = "prob")[,
    2]
hist(hbwip_trainpred, 50)

plot(density(hbwip_trainpred, na.rm = T, from = 0, to = 1))


terra::hist(hbwip, "WET", maxcell = 1e+07, breaks = 50)
## Warning: [hist] a sample of17% of the cells was used (of which 5% was NA)
terra::hist(hbwip_noriv, "WET", maxcell = 1e+07, breaks = 50)
## Warning: [hist] a sample of17% of the cells was used (of which 5% was NA)

hb_pedons <- vect(read.csv("UplandWetlandGradient/data/HBEF_pedon_locations.csv"),
    geom = c("easting", "northing"), crs = "EPSG:26919")

hb_pedons_ext <- as.data.frame(terra::extract(hbwip, hb_pedons,
    bind = T, method = "simple"))
hb_pedons_ext$hpu <- as.factor(hb_pedons_ext$hpu)

hb_pedons_ext |>
    group_by(hpu) |>
    reframe(result = mean(WET, na.rm = T))
## # A tibble: 9 × 2
##   hpu   result
##   <fct>  <dbl>
## 1 Bh    0.148 
## 2 Bhs   0.0447
## 3 Bhsm  0.0616
## 4 Bi    0.0750
## 5 E     0.0643
## 6 H     0.503 
## 7 I     0.193 
## 8 O     0.0280
## 9 T     0.0612

plot(hbwip$WET)
plot(hb_pedons[hb_pedons$hpu == "H"], "hpu", type = "classes",
    add = T)
plot(hbwip$WET)
plot(hb_pedons[hb_pedons$hpu == "Bh"], "hpu", type = "classes",
    add = T)

hb_pedons_ext |>
    arrange(WET) |>
    na.omit() |>
    mutate(hpu = factor(hpu, levels = c("H", "I", "Bh", "T",
        "Bi", "E", "Bhs", "Bhsm", "O"))) |>
    ggplot(aes(x = hpu, y = WET)) + geom_violin(draw_quantiles = c(0.5),
    adjust = 1.5) + geom_point(aes(color = hpu), alpha = 0.5,
    position = position_jitter(seed = 11))

ggplot(hb_pedons_ext, aes(x = reorder(pm, desc(WET), FUN = "median"),
    y = WET)) + geom_violin(draw_quantiles = c(0.25, 0.5, 0.75),
    adjust = 1.5) + geom_point(aes(color = pm), alpha = 0.5,
    position = position_jitter(seed = 11))
## Warning: Removed 4 rows containing non-finite values (`stat_ydensity()`).
## Warning: Groups with fewer than two data points have been dropped.
## Warning: Removed 4 rows containing missing values (`geom_point()`).

Is there a model?

hbmod <- lm(WET ~ hpu, data = hb_pedons_ext)
summary(hbmod)
## 
## Call:
## lm(formula = WET ~ hpu, data = hb_pedons_ext)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.19318 -0.06195 -0.04470  0.00085  0.74792 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.14808    0.01844   8.029 1.55e-14 ***
## hpuBhs      -0.10339    0.02941  -3.515 0.000498 ***
## hpuBhsm     -0.08648    0.05249  -1.648 0.100330    
## hpuBi       -0.07308    0.02750  -2.657 0.008246 ** 
## hpuE        -0.08379    0.03514  -2.385 0.017634 *  
## hpuH         0.35492    0.11142   3.185 0.001577 ** 
## hpuI         0.04509    0.04196   1.075 0.283294    
## hpuO        -0.12008    0.03934  -3.052 0.002445 ** 
## hpuT        -0.08693    0.02392  -3.634 0.000322 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1554 on 346 degrees of freedom
##   (4 observations deleted due to missingness)
## Multiple R-squared:  0.1119, Adjusted R-squared:  0.09132 
## F-statistic: 5.447 on 8 and 346 DF,  p-value: 1.759e-06
anova(hbmod)
## Analysis of Variance Table
## 
## Response: WET
##            Df Sum Sq Mean Sq F value    Pr(>F)    
## hpu         8 1.0524 0.13155  5.4471 1.759e-06 ***
## Residuals 346 8.3558 0.02415                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
w3_pedons <- vect(read.csv("UplandWetlandGradient/data/hbef_w3_lat_pedons/W3-lateral-weathering-pedons.csv"),
    geom = c("easting", "northing"), crs = "EPSG:26919")

w3_pedons_ext <- as.data.frame(terra::extract(hbwip, w3_pedons,
    bind = T, method = "bilinear"))

w3 <- hb[hb$WS == "WS3"]
w3_wip <- crop(hbwip, w3, mask = T)
plot(hbwip$WET, legend = F, ext = w3)
plot(w3_pedons, "hpu", type = "classes", add = T, legend = T,
    cex = 0.7)
lines(w3)
lines(hbnwi_hydrog)

ggplot(w3_pedons_ext, aes(x = reorder(hpu, desc(WET), FUN = "max"),
    y = WET)) + geom_violin(draw_quantiles = c(0.25, 0.5, 0.75),
    adjust = 1.5) + geom_point(aes(color = hpu), alpha = 0.5,
    position = position_jitter(seed = 11))
## Warning: Removed 1 rows containing non-finite values (`stat_ydensity()`).
## Warning: Removed 1 rows containing missing values (`geom_point()`).

ggplot(w3_pedons_ext, aes(x = water_table, y = WET)) + geom_point(aes(color = hpu),
    alpha = 0.5, position = position_jitter(seed = 11)) + stat_smooth(method = "lm",
    geom = "smooth", se = F)
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 30 rows containing non-finite values (`stat_smooth()`).
## Warning: Removed 30 rows containing missing values (`geom_point()`).

Is there a model in W3?

w3mod <- lm(WET ~ hpu, data = w3_pedons_ext)
summary(w3mod)
## 
## Call:
## lm(formula = WET ~ hpu, data = w3_pedons_ext)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.12298 -0.04011 -0.01673  0.01687  0.25955 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)   
## (Intercept)         0.116994   0.034248   3.416  0.00184 **
## hpuBhs podzol      -0.066107   0.048434  -1.365  0.18244   
## hpuBimodal podzol  -0.087680   0.076582  -1.145  0.26130   
## hpuE podzol        -0.078464   0.050134  -1.565  0.12805   
## hpuLithic Histosol -0.088856   0.052315  -1.698  0.09977 . 
## hpuTypical podzol   0.007971   0.055224   0.144  0.88620   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.09687 on 30 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.1668, Adjusted R-squared:  0.02797 
## F-statistic: 1.201 on 5 and 30 DF,  p-value: 0.3323
anova(w3mod)
## Analysis of Variance Table
## 
## Response: WET
##           Df  Sum Sq   Mean Sq F value Pr(>F)
## hpu        5 0.05637 0.0112739  1.2015 0.3323
## Residuals 30 0.28151 0.0093836
w3mod <- lm(WET ~ water_table, data = w3_pedons_ext)
summary(w3mod)
## 
## Call:
## lm(formula = WET ~ water_table, data = w3_pedons_ext)
## 
## Residuals:
##        3        6       11       17       19       28       35 
## -0.11031 -0.06748 -0.04230  0.02109  0.25246 -0.03985 -0.01360 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)
## (Intercept)  0.1561297  0.1427138   1.094    0.324
## water_table -0.0008764  0.0024205  -0.362    0.732
## 
## Residual standard error: 0.13 on 5 degrees of freedom
##   (30 observations deleted due to missingness)
## Multiple R-squared:  0.02555,    Adjusted R-squared:  -0.1693 
## F-statistic: 0.1311 on 1 and 5 DF,  p-value: 0.7321
anova(w3mod)
## Analysis of Variance Table
## 
## Response: WET
##             Df   Sum Sq   Mean Sq F value Pr(>F)
## water_table  1 0.002215 0.0022146  0.1311 0.7321
## Residuals    5 0.084466 0.0168931
hbmoistgrad <- read.csv("UplandWetlandGradient/data/hourly_soil_moisture_gradient_plots.csv")

plot(hbmoistgrad$E12_VWC_1)